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72 Abstract. We present an ab-initio study of neutron drops. We use Quantum Monte Carlo 

techniques to calculate the energy up to 54 neutrons in different external potentials, and we 
I compare the results with Skyrme forces. We also calculate the rms radii and radial densities, 

and we find that a re-adjustment of the gradient term in Skyrme is needed in order to reproduce 
the properties of these systems given by the ab-initio calculation. By using the ab-initio results 
for neutron drops for close- and open-shell configurations, we suggest how to improve Skyrme 

l"i forces when dealing with systems with large isospin-asymmetries like neutron-rich nuclei. 
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1. Introduction 

I Even though neutron drops, neutrons confined by an external field, cannot be studied 

experimentally, they have nevertheless generated intense theoretical interests for a couple of 
^vq reasons. Neutron drops provide a very simple model of neutron-rich nuclei, in which the core 

is modeled as an external potential acting on valence neutrons. For example in Refs. [H [2j |3] 
neutron-rich oxygen isotopes have been modeled by neutron drops, and in Ref. [I] the same model 
has been used to study calcium isotopes. Second, they provide a description of inhomogeneous 
neutron matter that can be used as data for calibrating model energy density functionals in 
several conditions El E] • 

On the experimental side, new facilities plan to study the physics of neutron-rich nuclei 
close to the drip line. The study of these systems is typically done by using methods based on 
energy density functionals and effective forces. These forces, like Skyrme or Gogny, are fitted 
to reproduce properties of nuclear and neutron matter combined with experimental energies 
of nuclei close to stability. By constraining the density functionals to the equation of state of 
nuclear and neutron matter, the symmetry energy is implicitly included in these forces. However, 
the terms dependent on derivatives of proton and neutron densities are constrained only to nuclei 
close to stability, and then for very small isospin-asymmetries. 

The use of these functionals to study nuclei close to the neutron drip line requires then an 
important extrapolation to large isospin-asymmetries. This extrapolation is even more drastic 
when the Skyrme forces are used to study the properties of the neutron star crust, where the 
matter is made by extremely neutron-rich nuclei surrounded by a sea of neutrons. For these 
reasons, ab-initio calculations of these systems starting from accurate nuclear Hamiltonians are 
important to constrain density functionals. 

Important advances have been made in building very accurate Hamiltonians describing 
nuclear systems. The modern nucleon-nucleon forces fit scattering data with high accuracy [El [9], 



and also the three-body forces are quite well constrained to reproduce properties of light 
nuclei jTO], [H] . 

Furthermore, in recent years new techniques have been developed to solve for the ground- 
state of nuclear systems. In this paper we present results obtained with Quantum Monte Carlo 
techniques, i.e. the well known Green's Function Monte Carlo (GFMC) and Auxiliary Field 
Diffusion Monte Carlo (AFDMC). The two different algorithms use different variational wave 
functions, but they show similar accuracy in calculating the energy of pure neutron systems. 



2. Nuclear Hamiltonian and Quantum Monte Carlo method 

In our model, neutrons are non-relativistic point-like particles interacting via two- and three- 
body forces, that are confined by an external potential: 

A 2 

# = Ef^ + E v -^) + E%' + E C 1 ) 

i=l i i<j i<j<k 

where the external potential is a Harmonic Oscillator (HO) 

V H o{r) = \muo 2 r\ (2) 

or a Woods-Saxon 

V ws {r) = -Vo— ^ r-^ . (3) 

1 + exp[(r - r )/aj 

In our study the parameters are Vo = —35.5 MeV, ro = 3 fm and a = 1.1 fm, and ft 2 /m = 41.44 
MeV-fm 2 . The parametrization of the WS has been chosen to reproduce the properties of oxygen 
isotopes pQ, and we consider different values of fiu to change the density of the system. 

The two body-potential that we use is the Argonne AV8' [12] , a simplified form of the Argonne 
AV18 [8j. Although simpler to use in QMC calculations, the AV8' provides almost the same 
accuracy as AV18 in fitting NN scattering data. 

In this work we use the UIX three-body force, that has been originally proposed in 
combination with the Argonne AV18 and AV8' |13| . Although it slightly underbinds the 
energy of light nuclei, it has been extensively used to study the equation of state of nuclear 
and neutron matter [TH [JS], and gives an equation of state of pure neutron matter stiff 
enough to support a neutron star with mass higher than 2 solar masses [HI [JSJ [THJ, [T7], as 
suggested by recent observations [18j [19] . The Illinois three-body forces have been introduced 
to improve the description of both ground- and excited-states of light nuclei with A < 8 |10j . 
and produce energy spectra for A up to 12 in excellent agreement with data [11], but produce 
an unphysical overbinding in pure neutron systems [20] . In neutron drops their contribution is 
also attractive [2T] . 

We solve the many-body ground-state using both the Green's Function Monte Carlo (GFMC) 
and Auxiliary Field Diffusion Monte Carlo (AFDMC) [22]. The main idea of QMC methods is 
to evolve a many-body wave function in imaginary-time: 

V(T)=exp[-HT}V v , (4) 

where *$> v is a variational ansatz, and H is the Hamiltonian of the system. In the limit of r — > oo, 
^ approaches the ground-state of H. The evolution in imaginary-time is performed by sampling 
configurations of the system using Monte Carlo techniques, and expectation values are evaluated 
over the sampled configurations. For more details see for example Refs. [15} 123]. 

The GFMC has proven to be extremely accurate to study properties of light nuclei. The 
variational wave function includes all the possible spin/isospin states of the nucleons and provides 
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Figure 1. The energy of neutrons in a HO well with %uj = 10 MeV (upper panel) and 5 MeV 
(lower panel) in units of ujN a ^ . The red dots are the results given by AFDMC, blue squares 
are from GFMC, and the black line is the result obtained using the Skyrme SLY4. The violet 
line is the adjusted SLY4 where the strength of the gradient, pairing, and spin-orbit terms have 
been changed. The figure is taken from Ref. [S]. 



a good variational ansatz to start the projection in imaginary time. The exponential growing 
of the spin-isospin states limits the calculation to the 12 C [2]. The AFDMC method does not 
include all the spin/isospin states in the wave function. The latter are instead sampled using 
the Hubbard-Stratonovich transformation. Using the AFDMC the calculation can be extended 
up to many neutrons, making the simulation of large neutron drops with up to several tens of 
neutrons possible. In the case of homogeneous systems the AFDMC has been used to simulate 
up to 114 neutrons [15], and up to 70 neutrons in an external well. 

3. Results 

The energy of neutron drops confined by Vho is shown in Fig. [I] for two different frequencies of 
the external potential. The red points are the results obtained using the AFDMC method, and 
the blue ones using the GFMC. The two solid lines are the results given by using the original 
Skyrme SLY4 force |24j . and a modified version. The energy is in units of the Thomas- Fermi 
energy, that is proportional to wY 4 / 3 , to see the extrapolation to the thermodynamic limit. The 
two QMC methods agree within 1% for the Tiuj = 10 MeV trap, and the difference increases 
up to 4% for Tiuj = 5 MeV. The reason is that pairing correlations are included in the GFMC 
wave function, but not in the AFDMC. At low densities neutrons are superfluids, and pairing 
correlations are quite important to include for open-shell configurations. 

The difference between QMC and Skyrme at closed shells is mainly due to two effects, the 
bulk contribution and the gradient term. Skyrme forces typically give an EOS of pure neutron 
matter at densities lower than saturation that is more repulsive than microscopic calculations. 
The equation of state of pure neutron matter is shown in Fig. [2] where we compare the AFDMC 
results from Refs. [13 [23 [29] , the GFMC calculation of Gezerlis and Carlson [261 EZ], and the 
equation of state given by SLY4. We make the reasonable assumption that Skyrme's bulk term 
cannot explain the difference between QMC and Skyrme energy in neutron drops. Then, since 
the pairing and the spin-orbit terms are expected to be very weak with respect to the gradient 
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Figure 2. The equation of state of pure neutron matter. The AFDMC results |15t I25j are 
compared to the QMC results at low densities of Gezerlis and Carlson |26t I27j . and to the 
equation of state of Skyrme SLY4 |28j . 



term for closed shell configurations, we can use the energy at N=8, 20 and 40 to re-adjust the 
gradient term of Skyrme. The energy of neutron drops with N near closed shells can be used 
to adjust the spin-orbit strength because for these configurations the pairing is not important. 
Finally, by comparing the energy of half-filled shells, we can tune the pairing term. 

The ab-initio results can also be used as a benchmark for new Skyrme forces. For example, 
in a recent paper Kortelainen et al. have derived new Skyrme forces, and they have produced 
neutron drops results as a prediction to be compared with QMC calculations. The comparison 
is shown in Fig. [3| The agreement between the two new Skyrme forces with QMC results is 
remarkable considering that the ab-initio energies of neutron drops have not been included as a 
constraint in producing the UNEDFO and UNEDF1 forces. 

The use of a Woods-Saxon external potential Vws is a ls° interesting because it reflects 
properties of nuclei. It saturates for some number of neutrons, providing a good test of the 
surface effects of Skyrme. The comparison between the original and the adjusted Skyrme SLY4 
with QMC for neutrons in a WS potential is shown in Fig. [4} We stress that the SLY4 has been 
adjusted only to reproduce the results in the HO well, and is used here to calculate the neutron 
drop energies as a prediction. Also in this case the agreement with the QMC calculations is 
quite good. 

Finally, we studied the effect of the adjusted Skyrme to the structure of neutron drops by 
calculating the rms radii and density. The rms radius V < r 2 > is shown in Fig. [5] where we 
compare the GFMC estimate with the original and adjusted SLY4 for different numbers of 
neutrons in the WS well (lower panel), and similarly with HO well for two different frequencies 
(upper panel). In Fig. [6] we show the radial density calculated with GFMC and Skyrme in the 
case of HO well for a closed shell configuration with N=8 (upper panel) , and for an open-shell 
with N=14 (lower panel). Both for rms radii and for the density the adjusted Skyrme agree 
much better with QMC calculations, showing that the energy constraint is sufficient to improve 
the accuracy of density functionals, whose predictive power for other operators would drastically 
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Figure 3. Energy of neutrons in a HO well. The AFDMC results are compared with those 
given by UNEDFO and UNEDF1 calculations. Figure taken from Ref. [30] with courtesy of W. 
Nazarewicz. 




Figure 4. The energy per particle of neutrons in a WS well obtained using AFDMC (blue 
squares) and GFMC (red points) compared to Skyrme. Figure taken from Ref. [5]. 
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Figure 5. The radii of neutron drops as a function of the number of neutrons calculated using 
GFMC (blue points) compared with the original (black line) and the adjusted SLY4 (violet line). 
In the top panel we show the results in the HO external potential, and in the bottom panel those 
obtained in a WS well. Figure taken from Ref. [5]. 



improve. 

4. Conclusions 

We performed ab-initio calculations of neutrons confined in external potentials. These systems 
are interesting to study because they provide a simple model to describe neutron-rich nuclei, and 
can be used to constrain energy density functionals. By using QMC techniques, we calculated 
the energy of neutrons in different potentials. We have considered both harmonic oscillator and 
Woods-Saxon wells in order to provide results for systems with rather different geometries. We 
paid particular attention to the comparison between the GFMC and AFDMC methods. The 
agreement is quite good, at the few percent level for the total energy of the system. 

We compared the energies of neutron drops given by QMC methods with those given by 
different Skyrme forces, and found that they always produce lower energies. Since the density 
of these systems is lower than nuclear densities, we suggest that the overbinding is mainly due 
to the (neutron) gradient term. We found that, in order to reproduce our ab-initio results, the 
pairing and the spin-orbit terms of Skyrme need to be changed too. 

We have used AFDMC to calculate the energy up to 54 neutrons, and GFMC to calculate the 
rms radii and radial densities of neutron drops. After re-adjustment of the SLY4 parametrization, 
the Skyrme force provides much better agreement for these quantities. 

Acknowledgments 

The author would like to thank J. Carlson and R. Schiavilla for critical comments on the 
manuscript, and W. Nazarewicz for the useful discussions and for the permission to show Fig.|3j 
This work is supported by DOE Grants No. DE-FC02-07ER41457 (UNEDF SciDAC) and No. 
DE-AC52-06NA25396, and by the LANL LDRD program. Computer time was made available by 
Los Alamos Open Supercomputing, and by the National Energy Research Scientific Computing 




r (fm) 



Figure 6. The radial density of neutron drops for different configurations calculated using 
GFMC (squares) compared with the original (dashed lines) and the adjusted SLY4 (solid lines) 
in a HO well. In the top panel we show the results of a closed shell configuration with N=8, and 
in the bottom panel the density of an open shell with N=14. The two colors are for different 
HO frequencies. Figure taken from Ref. [5]. 
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